Predicting the coherence resonance curve using a semi-analytical treatment 
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Emergence of noise induced regularity or Coherence Resonance in nonlinear excitable systems is well known. 
We explain theoretically why the normalized variance (Vn) of inter spike time intervals, which is a measure of 
regularity in such systems, has a unimodal profile. Our semi-analytic treatment of the associated spiking process 
produces a general yet simple formula for Vn, which we show is in very good agreement with numerics in two 
test cases, namely the FitzHugh-Nagumo model and the Chemical Oscillator model. 
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INTRODUCTION 

Many deterministic, nonlinear, excitable systems, for ex- 
ample, the FitzHugh-Nagumo model (FHN) or the Chem- 
ical Oscillator model (CO) 101, undergo bifurcation from a 
stable focus to a stable limit cycle (LC) behavior when a sys- 
tem parameter is tuned. However, holding the parameter near 
the bifurcation point, on the stable focus side the system can 
still be made to exhibit spiking behavior (which is otherwise 
the signature of a limit cycle), by adding a random uncorre- 
cted noise to the system. The noise forces the system to in- 
termittently jump across the bifurcation point in the parameter 
space. As a result of these random excursions, the system ex- 
hibits intermittent cyclic behavior which manifests as spikes 
in the dynamical variable. Interestingly, the time intervals Tp, 
between two successive noise driven spikes, which are in gen- 
eral irregular, strangely becomes fairly regular at an optimal 
noise value (the resonance point). This phenomenon is called 
Coherence Resonance. It has attracted considerable interest 
theoretically as well as experimentally lUflBElTllilgifloll. 
as quite counter-intuitively order arises with the aid of tuned 
randomness. A quantitative means of detecting this resonance 
point is enumerating the normalized variance (Vn) defined by 

Vn - 



'('''p) ~ {^p) I ("^p)' a function of noise strength. 
Here (.) denotes statistical time average. Typically Vn is enu- 
merated from time-series analysis of spikes generated by the 
system, subjected to noise. The noise strength at which mini- 
mum of Vn occurs is the desired point of resonance. 

The analytical work so far on this subject, have either dealt 
with a toy model Ull, or addressed special limits of the FHN 
model e.g, very weak noise 111 111 , and infinite time scale sepa- 
ration between the fast and slow variables S 12, 13]. A pio- 
neering qualitative understanding of the phenomenon is given 
by Pikovsky and Kurths ||l|], who argue that the resonance hap- 
pens as a competition between two time scales - the activation 
time ta (the time between the end of one spike and beginning 
of another) and the excursion time i.e., duration of a spike. 
The inter spike interval (ISI) Tp ^ ta + te- They claim that 
ta has a strong dependence on noise intensity and follows a 
simple Kramer's lll411 like formula, whereas te has a much 
weaker noise dependence and corresponds to the decay time 
of unstable excited state. Kramers theory describes the noise 



driven escape time Tcsc of a particle (say y) from a deep po- 
tential trap, and gives Tcsc ^ cxp{Eb/ D^); here D is noise 
amplitude, and Eb is the barrier height. But excitable systems 
with two coupled variables x and y pose new challenges: the 
barrier Eb is both dynamic and D dependent. The effective 
barrier for y is dynamic as it is generated by x which itself 
is a dynamical variable. Furthermore our numerical studies 
show that barrier parameters, like its width 6, are indeed D 
dependent. In this paper, we avoid invoking Kramers picture 
apriori, and show that the timescales and ta can be under- 
stood from alternative arguments. 

We derive below a simple theoretical formula for Vn, 
which will be generally applicable to any nonlinear system 
exhibiting coherence resonance. There are parameters in the 
universal formula, which depend on the specific details of 
the nonlinear system at hand, and can only be fixed by some 
amount of numerical or alternatively experimental analysis. 
Thus the formula is semi-theoretical. Although this may seem 
as no less work than the usual time-series analysis, as we show 
below, it certainly involves incorporation of enhanced under- 
standing of the phenomenon compared to what existed before. 
To support our claim of generality, we study two very differ- 
ent nonlinear systems: the FHN model HI] and the CO model 
iS [13 15]. We show that our predicted formula fits quite 
well, with the curve of Vn obtained by brute force time-series 
analysis, in both the cases. 



MODEL 

Before starting our main analysis, let us define the FHN and 
CO systems in the presence of noise, to make this paper self 
contained. The FHN model has the following equations 



dx 
'~dt 



y, 



dy 
dt 



x + a + Di{t). 



(1) 



Here a, D and e (<C 1) are the three parameters. For \a\ > 1, 
there is a stable fixed point at = —a , y^ = ^ — a, while 
for \a\ < la limit cycle exists in the x — y space and dy- 
namics of both the variables are periodic. The value of a 
on the fixed point side, which we hold fixed for our simu- 
lation, is denoted by oq. The parameter D is the amplitude 
of the Gaussian white noise ^, for which {£,{t)) = and 
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FIG. 1: Semi-log plot of Pi(rp) against Tp for FHN model. Here 
ao — 1.1 and e = 0.01 (see Eq. l[Tll) — these same values are 
used for other figures in the paper. The straight lines are exponential 
fits (placed higher for visual clarity) to the tail to obtain t"™ values 
and they are 2.225, 1.007 for D = 0.05 and 0.12 respectively. The 
corresponding (defined in the text) values are 3.573 and 2.904. 



(^(t)^(t )} = S{t — t ). The small parameter e makes the mo- 
tion on the limit cycle much faster along the x direction than 
the y. The second model of CO is defined by the following 
equations: 



V — u 



du 

''Itt ~ R 
dc u — V 

di " 



R 



- /(m,c), 

+ {l-c) + af{u,c) (2) 



where f{u, c) = c{aiu + a2U^ + a^u'^) and v = vq + D^{t). 
Here u and c are the dynamical variables and R, vq, ai, 02, 
03, e, a, and D are the parameters, v is the bifurcation param- 
eter. Limit cycle exists for the values v < 29.235 whereas 
for V > 29.235, a steady state fixed point behavior is ob- 
served. The system variables and parameters are derived from 
the reaction-rate kinetics of the interacting chemical species. 
The details regarding the construction of the model equation 
are furnished elsewhere 



RESULTS 

If one makes a simple-minded first guess that the inter spike 
intervals Tp have a Poisson distribution, then Vn would be a 
constant (independent of noise strength) which is empirically 
not the case. So what is the distribution of Tp? For a random 
train of spikes which are almost independent, it seems very 
likely that the distribution of ISl will have an exponential tail 
lfl6ll . Yet a specialty of the spikes in the non-linear systems of 
our concern, is that a new spike cannot arise until the last spike 
subsides. Thus Tp cannot be any smaller than characteristic 
'spike width' Wg (a finite quantity), i.e. the distribution of Tp 
is expected to have a sharp lower cutoff at some finite Tmin- 



We stress here that if this lower cutoff were absent, then Vn 
would have had no variation and coherence resonance would 
have vanished. Thus we expect the probability density of Tp 
to be. 



Pl(Tp) = NQ{Tp - T,r,iniD)) CXp 



(3) 



Here Q is the Heaviside function Ill7ll . while Tcsc is the char- 



acteristic time associated with exponential tail of Pi(Tp). In 
Eq. (O), the normalization constant N = (e'^""/'^''='= ) / Tcsc- We 
have checked that the distribution of Tp obtained from the time 
series analysis of the FHN and CO models are consistent with 
Eq. (O — see Fig. [T]for numerically obtained Pi (Tp) for the 
FHN system for two different D values. Despite the two val- 
ues of D, one being away and another close to the resonance 
point, one can see clearly that the shape of the curves Pi (Tp) 
shows no qualitative variation. Of course the quantities t"™ 
and T"g"™ (where the superscript "num" denotes numerical) 
are functions of D; in fact both decrease with D. The nota- 
tional distinction between t^™ in Fig. [T] and Tmin in Eq. ^ 
is necessary, as the numerical curve in Fig. [T] does not rise 
strictly as a 8 function. To be precise, in Fig.[T] t'^^",™ is de- 
fined as the average of the time t™,™ at which Pi(Tp) just 
starts becoming nonzero and the time Tj™\. at which Pi{Tp) 
reaches a peak. On the other hand r"g"™ is obtained by fitting 
an exponential to the tail of Pi(Tp). In this paper we attempt 
to obtain r,nin and Tosc theoretically, as opposed to the numer- 



and 



r"™ just described. Note that the 
" are analogous to the quantities Tg 



ical estimates r^^^^^ 
quantities t"™ and 
and Ta respectively as discussed in IJJ 

The first and the second moments of Tp, namely (rp) and 
(Tp), can be easily obtained using Eq. (O and using them in 
the definition of V^v we get 



Vn = 



,{D)+T,sciD)' 



(4) 



The simple formula for Vjy above, is the central result of this 
paper llSll . and is a good approximation in general for any 
non-linear system exhibiting coherence resonance, provided 
one could predict Tcsc{D) and T,„in(Z)) theoretically. In what 
follows we try to do the latter A similar formula as Eq. 

was derived, although in the low D limit ll ill for anti- 
coherence resonance. 

Formally, the resonance point is obtained by setting the 
derivative of Vjv w.r.t. D equals 0. That implies the following 
relation 



Here Dres denotes the value of D at the minimum of the Vn 
curve i.e. at the resonance point. T^^^{Dres) '^^'i '''nuni^ res) 
denote their respective derivatives with D evaluated at Dres 
. However, since both Tmin {Dres ) and Tcsc {Dres ) are system 
specific and are obtained partly numerically, the scope of the 
analytical application of Eq. ^ is limited. 
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FIG. 2; Top frame; rmin(I^) versus acff(D) in FHN. Bottom frame: 
-rmin(£') versus Weff(D) in CO. Here vo = 29.24, e = 0.03, a = 
0.1, ai = 1.125, a2 = -0.075, as = 0.00125, i? = 10 (see Eq. 
^) — these same values are used for CO in other figures. For both 
frames: The empty symbols are for t^^\ and filled symbols are for 
r^™i. The solid lines represent rmin from Eq. The numerical 
data and the theoretical curves shows excellent agreement. 



We start with a hypothesis about the functional dependence 
of Tmin on D. We claim that the action of noise on Eq. ([T]) (or 
|2| merely shifts a (or v) to UcS (or Vcs), with Oefi = a^ — D (or 
■^eff = vq — D). To be brief let us focus on the FHN system and 
the parameter a. The parameter value ag corresponds to the 
initial stable fixed point. The Ocff makes the system feel that 
it is on the LC side, across the bifurcation threshold ath = 1, 
and lead to a spike. The width of the spike Tmin is expected 
to be equal to the time period of the effective LC experienced, 
say ticp, i.e.. 



Tmin{D) =tiep(acff), and analogously for w. 



(6) 



Here we assume that t\cp, which is the property of the sys- 
tem is known apriori as a function of a. Note that the sys- 
tem can spike even if Ocff does not cross ath (and Tj^"™ can 
be measured numerically), but our above claim is not valid 
as t\cp is undefined. In the later case, we would claim that 
T'min(^) = Ws, the spikc width. 

We proceed to test Eq. (|6]l in FHN and CO models. In both 
the top (for FHN) and bottom (for CO) frames of Fig. |2] the 
solid lines are as per Eq. (|6]l. Instead of plotting r^""™, for 
more clarity, we have plotted t^™\ in empty symbols and 
r^™ in filled symbols. The fact that Tmin falls in between 
''"rrdrTr ''m™ ^'^^ range of D Studied, and the agreement 
being excellent for two distinct systems FHN and CO (with 
distinct t\cp{a) and t\cp{v) functions), gives strong empirical 
support for the formula in Eq. (|6]l. 

Next, we turn to Tcsc in Eq. (|4|i. The dynamics of one of 
the variables in the non-linear system, for example y in FHN 
or c in CO, under finite noise strength D, can be viewed as a 
stochastic process around the stable fixed points or c», re- 
spectively. For subsequent discussion we focus on FHN, but 
the results apply generally to any non-linear system exhibit- 
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FIG. 3: Tcsc{D) versus D for FHN (top frame) and CO (bottom 
frame) denoted by line joining filled symbols (the values of fo{D) 
and S,n(D) used for this plot are discussed in the text and fig. |4] 
Open symbols represents the numerical values r^c™(-D), obtained 
from fits as in fig[T] 



ing coherence resonance. Most often the noise displaces y a 
little and then it relaxes back to the fixed point, in a typical 
excursion time tq. Occasionally however, if the excursion of 
the variable (e.g. Ay ^ y — y^ in FHN) falls below a certain 
threshold denoted by a typical —5„i (here 5m > 0), the system 
exhibits a cycle and y exhibits a spike. The latter amounts to 
absorption of Ay at the boundary —S„i, on its first passage. 

Specific system dependent details of the shape of the ef- 
fective trapping potential is necessary to analytically calcu- 
late the above mentioned typical first passage time Tcsc- Since 
our purpose is to remain as general as possible, we make a 
simplifying general assumption that after every random kick 
the relaxation is instantaneous. In effect this is equivalent 
to coarse-graining in time over units of the typical excursion 
time To (mentioned above and defined below). 

Thus every excursion Ay 7^ at every discrete time step, 
may be treated as independent, and merely follows the noise 
and therefore has the same (Gaussian) distribution as the 
noise. Then it immediately follows, that the probability Q{n) 
that the signal Ay does not go below —Sm for n successive 
time steps and does so in the (n + 1)*'' step is 



Qin) = [P>i-6mTP<{-S,n) 



(7) 



whereP>(y) = -^/J^e^l^^da;' and P< = 1-F>. Eq© 
shows that Q{n) is exponential distributed, and its decay con- 
stant gives the "typical first passage time" lfl6ll in units of tq: 



P>{-Srn) = {l + erf{Sm/D))/2 



where 



(8) 



and erf (.) is the Error function jV, 

Note that the D dependence of Tosc comes from explicit 
dependence of P> on D, as well as the implicit dependence 
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FIG. 4: (Color online) (a) PD of to for the FHN along with exponential fit e"*''^" for two different D values, (b) PD of <5, and depiction of S,n 
as most probable 5 for the same two D as in (a), (c) fo versus D for FHN (o) and CO (A) extracted from figure (a). The dashed line segments 
are explained in the text, (d) Sm versus D extracted from figure (b) - two sets are for FHN (o) and CO (A). 



of the time unit tq and barrier location S„i on D. Of course 
To and 6m will be system specific and incorporate the detail 
nature of the dynamic potential trap. The procedure to find 
fo{D) and 6m{D) will be discussed later. If we assume that 
the latter two quantities are known a priori then Eq. (O maybe 
claimed to be a "theoretical" formula, and compared to the 
numerical values of r™™ obtained as in Fig. [T] In Fig. |3]we 
see that the agreement between the theoretical formula and 
numerical data are excellent. 

Using the asymptotic expansion of crf(.) IitIi in Eq. ^ 
we get Tcsc/tq ~ constant for D/6m ^ 1 and w e^^l^ 
for D /5„i <^ 1. The latter behavior has been referred to as 
Kramer's formula for Tosc id [HI, but one needs to be careful 
— unlike the usual Kramer's escape time formula, is not 
the barrier height of the potential well but rather proportional 
to the width of the well. 

What remains to be discussed is determination of tq{D) 
and Sm{D). To define tq precisely, we note that between 
two successive spikes of y, the process Ay (and Ac for CO 
model) crosses zero several times. Let tq be the time interval 
between zero crossings of Ay which is same as the excur- 
sion time mentioned earlier. A probability distribution (PD) 



of To is then found for every D, and the PD has an exponen- 
tial tail as shown in Fig. SJa). We define the time constant 
of the latter exponential fit to be fo{D). For FHN and CO 
systems the fo{D) thus obtained are shown in Fig.Hfc). But 
with increasing D the time stretches between two spikes be- 
come very small, making determination of fo(-D) unreliable 
due to poor statistics. So we took fQ{D) to be a constant (de- 
noted by dashed line segments in Fig. SJc)), for the D values 
beyond which fo{D) could not be reliably determined. A poi- 
renor; justification of the latter adhoc assumption for fo lies in 
the successful agreement with numerical data of Tcsc{D) (see 
Fig.EJ. 

Next, we define —S as the threshold of Ay at which spiking 
occurs. Then the PD of S for every D can be computed (see 
Fig. |4jb)) and the most probable value may be identified as 
6m- Plot of 6m is shown against D for both FHN and CO 
systems in Fig. HJd). These values of Sm were used to obtain 
the theoretical curve in Fig. [3] 

Finally, one can directly plot the Vat from the theoretical 
formulas in Eqs. (|4]i, (|6]l, and ^ and compare it with numer- 
ical Vn obtained from time series analysis (see Fig.|5]l. Both 
for FHN and CO the agreement is quite good and the locations 
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FIG. 5: Normalized Variance Plot for FHN(top) and CO(bottom). 
Open symbol denote the Vn from numerical time series analysis. 
While, filled symbols joined by line represents the theoretical values 
from Eq. The semi-analytic treatment was found to be valid for 
the following range of D values: For the FHN model 0.04 < D < 
0.2. For the CO model 0.05 < D < 0.6. 

of resonance (the minima) are obtained within acceptable er- 
ror limits. 

CONCLUSION 

Thus we claim to have found an alternate way of deter- 
mining Vn for non-linear systems exhibiting coherence res- 
onance, based on theoretical considerations rather than brute 
force time series analysis. Only three empirical inputs are re- 
quired for a specific system, namely (i) the Umit cycle period 
t\cp{a) as a function of the control parameter a, (ii) the typi- 
cal zero-crossing time interval 7^(Z?) of the relevant stochastic 
dynamical variable, and (iii) the typical distance of excursion 
6jn{D) beyond which the variable maybe regarded as "ab- 
sorbed" (i.e. it spikes). We highlight the fact that the effective 
barrier parameters fg and S„i turn out to be D dependent. It 



may seem no less work to obtain the empirical inputs (i)-(iii) 
for a system, yet once obtained they can be substituted in the 
simple theoretical formulas Eqs. (HJi, (|6]l, and ([8]l and coher- 
ence resonance maybe predicted. 
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